#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

extern double compute_curvature(),compute_delta(),
              compute_delta_general(),compute_heaviside(),
              compute_heaviside_general(),compute_delta_prime(),
              compute_heaviside_nozero(),compute_delta_nozero();
extern void get_current_area(),solid_body_flow();

#define mass_spread (density_spread)

#define maxeno 9
#define enodim 18
#define idx1 7
#define idx2 8
#define uplimit (enodim-maxeno)

#define all_non_conserve FALSE
#define shu_osher_flag FALSE
#define harten_flag FALSE

extern double b00[imax][jmax],b01[imax][jmax],b02[imax][jmax],
         b10[imax][jmax],b11[imax][jmax],b12[imax][jmax],
         b20[imax][jmax],b21[imax][jmax],b22[imax][jmax];
extern int number_time_steps;

static double ubar,vbar,wave_speed,uleft,uright,vbot,vtop;
static int ix,iy;
static logical xdir;
static double aa[enodim],bb[enodim],cc[enodim],dd[enodim],
              ee[enodim],ff[enodim],gg[enodim],hh[enodim];

static int eno_limit;

double table1(iv,inc,inc1)
int iv,inc,inc1;
{
int xindex,yindex;
double r_val;

  xindex=ix;
  yindex=iy;
  if (xdir)
   { xindex +=inc;yindex +=inc1; }
  else
   { yindex +=inc;xindex +=inc1; }

  r_val=1.0;
  if (use_axisymmetric) r_val=x[xindex];

  if ((xindex>=1)&&(xindex<=mx)&&(yindex>=1)&&(yindex<=my))
   {
    if (iv==1) return u[xindex][yindex];
    else if (iv==2) return v[xindex][yindex];
    else if (iv==3) return d1[xindex][yindex];
    else if (iv==4) return d2[xindex][yindex];
    else if (iv==5) return d3[xindex][yindex];
    else if (iv==6) 
      return r_val*u[xindex][yindex]*b20[xindex][yindex];
    else if (iv==7) return v[xindex][yindex]*b20[xindex][yindex];
    else if (iv==8)
     return r_val*u[xindex][yindex]*d1[xindex][yindex];
    else if (iv==9)
     return r_val*v[xindex][yindex]*d1[xindex][yindex];
    else if (iv==10)
     return r_val*u[xindex][yindex];
    else if (iv==11)
     return r_val*v[xindex][yindex];
   }
  else
   {
    if (iv==1) return gu(xindex,yindex);
    else if (iv==2) return gv(xindex,yindex);
    else if (iv==3) return gd11(xindex,yindex);
    else if (iv==4) return gd2(xindex,yindex);
    else if (iv==5) return gd3(xindex,yindex);
    else if (iv==6)
     return compute_heaviside_nozero(gd11(xindex,yindex),mass_spread)*
          gru(xindex,yindex);
    else if (iv==7)
     return compute_heaviside_nozero(gd11(xindex,yindex),mass_spread)*
          gv(xindex,yindex);
    else if (iv==8)
     return gru(xindex,yindex)*gd11(xindex,yindex);
    else if (iv==9)
     return grv(xindex,yindex)*gd11(xindex,yindex);
    else if (iv==10)
     return gru(xindex,yindex);
    else if (iv==11)
     return grv(xindex,yindex);
   }

}

double table(iv,inc)
int iv,inc;
{
  return table1(iv,inc,0);
}

double eno_conserve(iv,k0)
int iv,k0;
{
double c0,c1,c2,c3,c4,c5,c6,c7,c8,
       q0,q1,q2,q3,q4,q5,q6,q7,q8;
int k1,k2,k3,k4,k5,k6,k7,k8,i;

/* compute u_{i+1/2} depending on wind type k0 */

 for (i=maxeno-eno_limit;i<uplimit+eno_limit;i++)
  aa[i]=table(iv,i-idx1)-table(iv,i-idx2);
 for (i=maxeno-eno_limit;i<uplimit+eno_limit-1;i++)
  bb[i]=(aa[i+1]-aa[i])*0.5;
if (eno_limit>3)
{
 for (i=maxeno-eno_limit;i<uplimit+eno_limit-2;i++)
  cc[i]=(bb[i+1]-bb[i])/3.0;
if (eno_limit>4)
{
 for (i=maxeno-eno_limit;i<uplimit+eno_limit-3;i++)
  dd[i]=(cc[i+1]-cc[i])/4.0;
 for (i=maxeno-eno_limit;i<uplimit+eno_limit-4;i++)
  ee[i]=(dd[i+1]-dd[i])/5.0;
 for (i=maxeno-eno_limit;i<uplimit+eno_limit-5;i++)
  ff[i]=(ee[i+1]-ee[i])/6.0;
 for (i=maxeno-eno_limit;i<uplimit+eno_limit-6;i++)
  gg[i]=(ff[i+1]-ff[i])/7.0;
 for (i=maxeno-eno_limit;i<uplimit+eno_limit-7;i++)
  hh[i]=(gg[i+1]-gg[i])/8.0;
}    /* no need for big table if eno_limit <=4 */
}   /* no need for big table if eno_limit <=3 */

 q2=q3=q4=q5=q6=q7=q8=0.0;
 c2=c3=c4=c5=c6=c7=c8=0.0;
 c0=table(iv,k0);

 q1=(1-2*k0);

 if (fabs(aa[idx2+k0])<fabs(aa[idx1+k0]))
  {
   k1=k0;c1=aa[idx2+k0]/2;q2=q1*(-k0-1);
  }
 else
  {
   k1=k0-1;c1=aa[idx1+k0]/2;q2=q1*(-k0+2);
  }
if (eno_limit>=3)
{
 if (fabs(bb[k1+idx2])<fabs(bb[k1+idx1]))
  {
   k2=k1;c2=bb[k1+idx2]/3;q3=q2*(-k1-2);
  }
 else
  {
   k2=k1-1;c2=bb[k1+idx1]/3;q3=q2*(-k1+2);
  }

if (eno_limit>=4)
{
 if (fabs(cc[k2+idx2])<fabs(cc[k2+idx1]))
  {
   k3=k2;c3=cc[k2+idx2]/4;q4=q3*(-k2-3);
  }
 else
  {
   k3=k2-1;c3=cc[k2+idx1]/4;q4=q3*(-k2+2);
  }

if (eno_limit>=5)
{
 if (fabs(dd[k3+idx2])<fabs(dd[k3+idx1]))
  {
   k4=k3;c4=dd[k3+idx2]/5;q5=q4*(-k3-4);
  }
 else
  {
   k4=k3-1;c4=dd[k3+idx1]/5;q5=q4*(-k3+2);
  }

if (eno_limit>=6)
{
 if (fabs(ee[k4+idx2])<fabs(ee[k4+idx1]))
  {
   k5=k4;c5=ee[k4+idx2]/6;q6=q5*(-k4-5);
  }
 else
  {
   k5=k4-1;c5=ee[k4+idx1]/6;q6=q5*(-k4+2);
  }

if (eno_limit>=7)
{
 if (fabs(ff[k5+idx2])<fabs(ff[k5+idx1]))
  {
   k6=k5;c6=ff[k5+idx2]/7;q7=q6*(-k5-6);
  }
 else
  {
   k6=k5-1;c6=ff[k5+idx1]/7;q7=q6*(-k5+2);
  }

if (eno_limit>=8)
{
 if (fabs(gg[k6+idx2])<fabs(gg[k6+idx1]))
  {
   k7=k6;c7=gg[k6+idx2]/8;q8=q7*(-k6-7);
  }
 else
  {
   k7=k6-1;c7=gg[k6+idx1]/8;q8=q7*(-k6+2);
  }

if (eno_limit>=9)
{
 if (fabs(hh[k7+idx2])<fabs(hh[k7+idx1]))
  {
   k8=k7;c8=hh[k7+idx2]/9;
  }
 else
  {
   k8=k7-1;c8=hh[k7+idx1]/9;
  }
} /* 9 */ } /* 8 */ } /* 7 */ } /* 6 */ } /* 5 */
} /* 4 */ } /* 3 */

 q0=c0+c1*q1+c2*q2+c3*q3+c4*q4+c5*q5+c6*q6+c7*q7+c8*q8;
 return q0;
}   /* eno_conserve */

double eno_nonconserve(iv,i_upwind)
int iv,i_upwind;
{
double q0;
int k0,k1,k2,k3,k4,k5,k6,k7,k8,i,
  k1hold,k2hold,k3hold,k4hold;

/* compute u_x{i} depending on wind type k0 */

 for (i=0;i<10;i++)
  aa[i]=(table(iv,i-4)-table(iv,i-5))/h;
 for (i=0;i<9;i++)
  bb[i]=(aa[i+1]-aa[i])*0.5;
 for (i=0;i<8;i++)
  cc[i]=(bb[i+1]-bb[i])/3.0;
 for (i=0;i<7;i++)
  dd[i]=(cc[i+1]-cc[i])/4.0;
 for (i=0;i<6;i++)
  ee[i]=(dd[i+1]-dd[i])/5.0;

 k0=0;   /* non-conserve always include current point */
 k1=( (i_upwind==0) ? k0-1 : k0 );
 k1hold=-2*k1-1;

 if (fabs(bb[k1+4])<fabs(bb[k1+5]))
  {
   k2=k1-1;k2hold=-k2;
  }
 else
  {
   k2=k1;k2hold=-(k2+2);
  }
 if (fabs(cc[k2+4])<fabs(cc[k2+5]))
  {
   k3=k2-1;k3hold=-k3;
  }
 else
  {
   k3=k2;k3hold=-(k3+3);
  }
 if (fabs(dd[k3+4])<fabs(dd[k3+5]))
  {
   k4=k3-1;k4hold=-k4;
  }
 else
  {
   k4=k3;k4hold=-(k4+4);
  }
 if (fabs(ee[k4+4])<fabs(ee[k4+5]))
  k5=k4-1;
 else
  k5=k4;

 q0=aa[k1+5];
 if (eno_limit>1)
 q0+=bb[k2+5]*(k1hold); 
 if (eno_limit>2)
  q0+=cc[k3+5]*(k2hold)*(k1hold);
 if (eno_limit>3)
  q0+=dd[k4+5]*(k3hold)*(k2hold)*(k1hold);
 if (eno_limit>4)
  q0+=ee[k5+5]*(k4hold)*(k3hold)*(k2hold)*(k1hold);

 return q0;
}   /* eno_nonconserve */

double eno_edgeval(iv,i_upwind)
int iv,i_upwind;
{
double q0,coeffhold;
int k0,k1,k2,k3,k4,k5,k6,k7,k8,i,
  k0hold,k1hold,k2hold,k3hold,k4hold;

/* compute u_i+1/2 depending on wind type i_upwind */

 for (i=0;i<10;i++)
  aa[i]=(table(iv,i-4)-table(iv,i-5));
 for (i=0;i<9;i++)
  bb[i]=(aa[i+1]-aa[i])*0.5;
 for (i=0;i<8;i++)
  cc[i]=(bb[i+1]-bb[i])/3.0;
 for (i=0;i<7;i++)
  dd[i]=(cc[i+1]-cc[i])/4.0;
 for (i=0;i<6;i++)
  ee[i]=(dd[i+1]-dd[i])/5.0;

 k0=( (i_upwind==0) ? 0 : 1);  /* we want actual value */
 k0hold=-k0;

 if (fabs(aa[k0+4])<fabs(aa[k0+5]))
  { 
   k1=k0-1;k1hold=-k1;
  }
 else
  {
   k1=k0;k1hold=-(k1+1);
  }
 if (fabs(bb[k1+4])<fabs(bb[k1+5]))
  {
   k2=k1-1;k2hold=-k2;
  }
 else
  {
   k2=k1;k2hold=-(k2+2);
  }
 if (fabs(cc[k2+4])<fabs(cc[k2+5]))
  {
   k3=k2-1;k3hold=-k3;
  }
 else
  {
   k3=k2;k3hold=-(k3+3);
  }
 if (fabs(dd[k3+4])<fabs(dd[k3+5]))
  {
   k4=k3-1;k4hold=-k4;
  }
 else
  {
   k4=k3;k4hold=-(k4+4);
  }
 if (fabs(ee[k4+4])<fabs(ee[k4+5]))
  k5=k4-1;
 else
  k5=k4;

 q0=table(iv,k0);
 coeffhold=0.5+k0hold;
 q0+=aa[k1+5]*coeffhold;
 coeffhold*=(0.5+k1hold);
 q0+=bb[k2+5]*coeffhold;
 coeffhold*=(0.5+k2hold);
 if (eno_limit>2)
  q0+=cc[k3+5]*coeffhold;
 coeffhold*=(0.5+k3hold);
 if (eno_limit>3)
  q0+=dd[k4+5]*coeffhold;
 coeffhold*=(0.5+k4hold);
 if (eno_limit>4)
  q0+=ee[k5+5]*coeffhold;

 return q0;
}    /* eno_edgeval */


double enogeneral(iv,k0)
int iv,k0;
{
  if (all_non_conserve)
   return eno_nonconserve(iv,k0);
  else if (harten_flag)
   return eno_edgeval(iv,k0);
  else
  return eno_conserve(iv,k0);
}  /* enogeneral */

double enohalf1(iv)
int iv;
{
  return ( (wave_speed<=0.0) ? enogeneral(iv,1) : enogeneral(iv,0));
}

double eno_deriv(xhalf,yhalf)
double xhalf[imax][jmax],yhalf[imax][jmax];
{
double hold;

  if (all_non_conserve)
   return -ubar*xhalf[ix][iy]-vbar*yhalf[ix][iy];
  else if (harten_flag)
   {
    hold=(22.0*xhalf[ix][iy]+xhalf[ix][iy-1]+xhalf[ix][iy+1])-
         (22.0*xhalf[ix-1][iy]+xhalf[ix-1][iy-1]+xhalf[ix-1][iy+1])+
         (22.0*yhalf[ix][iy]+yhalf[ix-1][iy]+yhalf[ix+1][iy])-
         (22.0*yhalf[ix][iy-1]+yhalf[ix-1][iy-1]+yhalf[ix+1][iy-1]);
    return -(hold/(24.0*h));
   }
  else if (shu_osher_flag)
   return -(xhalf[ix][iy]-xhalf[ix-1][iy])/h-
           (yhalf[ix][iy]-yhalf[ix][iy-1])/h;
  else
   return -ubar*(xhalf[ix][iy]-xhalf[ix-1][iy])/h-
           vbar*(yhalf[ix][iy]-yhalf[ix][iy-1])/h;
}


/* compute convetion terms using "conservative" eno differencing.
   input ix_in,iy_in
   output cfl_ldiff - velocity (wave speed)
          ut,vt,d1,d2t,d3t
*/
void ldiff(cfl_ix,cfl_iy,cfl_ldiff)
double *cfl_ldiff;
int *cfl_ix,*cfl_iy;
{
double r_val,r_plus,r_minus,temp;
double uctrx,ubot,utop,uctry;
double vleft,vright,vctrx,vctry;
double d1left,d1right,d1ctrx,d1bot,d1top,d1ctry;
double d2left,d2right,d2ctrx,d2bot,d2top,d2ctry;
double d3left,d3right,d3ctrx,d3bot,d3top,d3ctry;
double area;
logical use_stream_fn;

if (use_debug)
printf("diff_type,eno_order,allnoncons.,so_flag,harten_flag %5d%5d%5d%5d%5d\n",
  diff_type,eno_order,all_non_conserve,shu_osher_flag,harten_flag);

  if (solid_body() )
   solid_body_flow();  /* fixed stream function/velocity field */
  else 
   for_all_points
    {
     if (! point1_ok)
      s[ix][iy]=gst1(s,ix,iy);
     if (! point_ok)
      {
       u[ix][iy]=gu(ix,iy);
       v[ix][iy]=gv(ix,iy);
      }
    }

  if (kill_vof_flag)
   for_domain
    b20[ix][iy]=compute_heaviside_nozero(d1[ix][iy],mass_spread);

  /* ix and iy are global to this module */
  for (iy=0;iy<=my;iy++)
   for (ix=0;ix<=mx;ix++)
    {
     r_plus=r_minus=r_val=h;
     if (use_axisymmetric)
      {
       r_val*=x[ix];
       r_plus*=(x[ix]+h/2.0);
      }

     xdir=TRUE;

     eno_limit=eno_order;

     if (shu_osher_flag)
      wave_speed=(u[ix][iy]+u[ix+1][iy])/2.0;
     else if (all_non_conserve)
      wave_speed=u[ix][iy];
     else if ((use_axisymmetric)&&(ix==0))
      wave_speed=0.0;
     else
      wave_speed=(u[ix][iy]+u[ix+1][iy])/2.0;
/*      wave_speed=(s[ix][iy]-s[ix][iy-1])/r_plus;*/

     if (! solid_body())
      {
       b00[ix][iy]=enohalf1(1);
       b02[ix][iy]=enohalf1(2);
      }

     if (harten_flag)
      b11[ix][iy]=enohalf1(3)*enohalf1(10);   
     else if (shu_osher_flag)
      b11[ix][iy]=enohalf1(8);
     else
      b11[ix][iy]=enohalf1(3);

     if (kill_vof_flag)
      b21[ix][iy]=enohalf1(6);

     if (three_phase() )
     {
      b20[ix][iy]=enohalf1(4);
      b22[ix][iy]=enohalf1(5);
     }

     xdir=FALSE;

     eno_limit=eno_order;

     if (shu_osher_flag)
      wave_speed=(v[ix][iy]+v[ix][iy+1])/2.0;
     else if (all_non_conserve)
      wave_speed=v[ix][iy];
     else
      wave_speed=(v[ix][iy]+v[ix][iy+1])/2.0;
/*      wave_speed=-(s[ix][iy]-s[ix-1][iy])/r_val; */

     if (! solid_body())
      {
       b10[ix][iy]=enohalf1(2);
       b01[ix][iy]=enohalf1(1);
      }

     if (harten_flag)
      b12[ix][iy]=enohalf1(3)*enohalf1(11);
     else if (shu_osher_flag)
      b12[ix][iy]=enohalf1(9);
     else
      b12[ix][iy]=enohalf1(3);

     if (kill_vof_flag)
      b22[ix][iy]=enohalf1(7);

     if (three_phase() )
     {
      b21[ix][iy]=enohalf1(4);
      d4[ix][iy]=enohalf1(5);
     }
    }    /* computing all the fluxes */ 

  *cfl_ldiff=0.0;
  *cfl_ix=*cfl_iy=1;

  for_domain
    {
     ubar=u[ix][iy];
     vbar=v[ix][iy];

     if (solid_body() )
      {
       ut[ix][iy]=0.0;
       vt[ix][iy]=0.0;
      }
     else
      {
       ut[ix][iy]=eno_deriv(b00,b01);
       vt[ix][iy]=eno_deriv(b02,b10);
      }

  
     d1t[ix][iy]=eno_deriv(b11,b12);
     if (exact_flag)  /* create d1t using a lagrangian formulation+interp */
      {
       d2t[ix][iy]=ubar;  /* will use later -- really is xt */
       d3t[ix][iy]=vbar;  /* yt */
      }

     if (three_phase() )
      {
       d2t[ix][iy]=eno_deriv(b20,b21);
       d3t[ix][iy]=eno_deriv(b22,d4);
      }

/* we compute change of volume for a PARTICULAR CELL,
   since uleft as assigned above is /r_val,
   we have (ur)_r = r(u)_r....   
   REMEMBER: r_val= h or rh 
*/

     if (kill_vof_flag)
      {
/* if doing integral style don't forget r_val */
       d5[ix][iy]=-r_val*( (b21[ix][iy]-b21[ix-1][iy])*h/r_val+
           (b22[ix][iy]-b22[ix][iy-1]) );
      }
      
     r_val=x[ix];
     if ((use_axisymmetric) && (r_val<1.0))
      ubar /= r_val;  
     temp=sqrt(sqr(ubar)+sqr(vbar));

     if (temp>*cfl_ldiff)
      {
       *cfl_ldiff=temp;
       *cfl_ix=ix;*cfl_iy=iy;
      }
    }    /* computing using precomputed fluxes */
}   /* ldiff */

void correct_advection(dd,ddt,ddt_hld,dd_hld)
double dd[imax][jmax],ddt[imax][jmax],
       ddt_hld[imax][jmax],dd_hld[imax][jmax];
{
double area,f_center,r_val,r_plus,r_minus,delta_prime,
  cl_int[3][3],bl_int[3][3],hnp[3][3],hn[3][3],
  bl_sum,cl_sum,lambda_linear,linear_error,quadratic_error,
  a_int[3][3],b_int[3][3],c_int[3][3],a_sum,b_sum,c_sum,delta_temp,
  disc,f_temp,lambda,lambda_1,lambda_2,non_ht,
  hnp_linear[3][3],hnp_quadratic[3][3],ht_linear,ht_quadratic;
int ix,iy,icx,icy,area_iterate;

#define max_area_it 1
/*
we just have phi_t=-H_x dx/H'

#define sum_int(var) \
 ( (r_minus*(var[0][0]+var[0][2])+r_plus*(var[2][2]+var[2][0])+ \
    6*(r_minus*var[0][1]+r_plus*var[2][1]+r_val*(var[1][0]+var[1][2]))+ \
    36*r_val*var[1][1])*sqr(h)/64.0 )
*/
#define sum_int(var) (var[1][1]*r_val)

  if (kill_vof_flag)
   {
    printf("correcting mass ...\n");
    for_domain
     {
      r_val=1.0;
      if (use_axisymmetric)
       r_val=x[ix];
      ddt[ix][iy]=d5[ix][iy]/
        (sqr(h)*r_val*compute_delta_nozero(dd[ix][iy],mass_spread));
     }   /* looping the domain */

   }   /* if kill_vof_flag */
}  /*   correct_advection */

